﻿using System;
using System.Collections.Generic;

namespace Allegro.MathLib.Functions
{
    public class Bessel
    {
        const int NUSE1=7, NUSE2=8;
        const int MAXIT = 50000;
        const double XMIN = 2.0;

	    double[] c1 , c2;
	    double xo,nuo;
	    double jo,yo,jpo,ypo;
	    //double io,ko,ipo,kpo;
	    //double aio,bio,aipo,bipo;
	    double sphjo,sphyo,sphjpo,sphypo;
	    int sphno;

        public double EPS = 1.0e-16;
        public double FPMIN;

	    public Bessel()
	    {
            //FPMIN = double.Epsilon / EPS;
	        FPMIN = 1.0e-30;

            xo = 9.99e99;
            nuo = 9.99e99;
            sphno = -9999;
	        c1 = new double[]
            {
                -1.142022680371168e0,
                6.5165112670737e-3,
	            3.087090173086e-4,
                -3.4706269649e-6,
                6.9437664e-9,
                3.67795e-11,
	            -1.356e-13
            };
            c2 = new double[]
            {
                1.843740587300905e0,
                -7.68528408447867e-2,
	            1.2719271366546e-3,
                -4.9717367042e-6,
                -3.31261198e-8,
                2.423096e-10,
	            -1.702e-13,
                -1.49e-15
            };
	    }

        /// <summary>
        /// Spherical Bessel function of order n
        /// </summary>
        /// <param name="n"></param>
        /// <param name="x"></param>
        /// <returns></returns>
        public double jl(int n, double x)
        {
            if (n != sphno || x != xo)
            {
                if (x > 0)
                {
                    sphbes(n, x);
                }
                else if (x == 0.0)
                {
                    if (n == 0)
                        sphjo = 1.0;
                    else
                        sphjo = 0.0;
                }
                else 
                    throw new ArgumentException("Invalid argument to Bessel.jl: x must be >= 0");
            }
            // The following is a dirty and very pragmatic way to handle the inability of the implementation to evaluate 
            // Bessel functions very close to zero. The problems arise when numbers in the besseljy method exceed 
            // max. or min. values for double precision numbers. 
            if (double.IsNaN(sphjo) && x < (double)n)
                sphjo = 0.0;
            return sphjo;
        }
        double sphbesy(int n, double x)
        {
            if (n != sphno || x != xo) sphbes(n, x);
            return sphyo;
        }

	    void besseljy(double nu, double x)
        {
	        double PI = System.Math.PI;
	        double a,b,br,bi,c,cr,ci,d,del,del1,den,di,dlr,dli,dr,e,f,fact,fact2,
		        fact3,ff,gam,gam1,gam2,gammi,gampl,h,p,pimu,pimu2,q,r,rjl,
		        rjl1,rjmu,rjp1,rjpl,rjtemp,ry1,rymu,rymup,rytemp,sum,sum1,
		        temp,w,x2,xi,xi2,xmu,xmu2,xx;
	        int i,isign,l,nl;

	        if (x <= 0.0 || nu < 0.0) 
                throw new ApplicationException("bad arguments in besseljy");

	        nl=(x < XMIN) ? (int) (nu+0.5) : (int)System.Math.Max(0,nu-x+1.5);
	        xmu=nu-nl;
	        xmu2=xmu*xmu;
	        xi=1.0/x;
	        xi2=2.0*xi;
	        w=xi2/PI;
	        isign=1;
	        h=nu*xi;
	        if (h < FPMIN) 
                h=FPMIN;
	        b=xi2*nu;
	        d=0.0;
	        c=h;
	        for (i=0; i<MAXIT; i++) 
            {
		        b += xi2;
		        d=b-d;
		        if (System.Math.Abs(d) < FPMIN) 
                    d=FPMIN;
		        c=b-1.0/c;
                if (System.Math.Abs(c) < FPMIN) 
                    c = FPMIN;
		        d=1.0/d;
		        del=c*d;
		        h=del*h;
		        if (d < 0.0) 
                    isign = -isign;
                if (System.Math.Abs(del - 1.0) <= EPS) 
                    break;
	        }
	        if (i >= MAXIT)
		        throw new ApplicationException("x too large in besseljy; try asymptotic expansion");
	        rjl=isign*FPMIN;
	        rjpl=h*rjl;
	        rjl1=rjl;
	        rjp1=rjpl;
	        fact=nu*xi;
	        for (l=nl-1; l>=0; l--) 
            {
		        rjtemp=fact*rjl+rjpl;
		        fact -= xi;
		        rjpl=fact*rjtemp-rjl;
		        rjl=rjtemp;
	        }
	        if (rjl == 0.0) 
                rjl=EPS;
	        f=rjpl/rjl;
	        if (x < XMIN) 
            {
		        x2=0.5*x;
		        pimu=PI*xmu;
                fact = (System.Math.Abs(pimu) < EPS ? 1.0 : pimu / System.Math.Sin(pimu));
		        d = -System.Math.Log(x2);
		        e=xmu*d;
                fact2 = (System.Math.Abs(e) < EPS ? 1.0 : System.Math.Sinh(e) / e);
		        xx=8.0*(xmu*xmu)-1.0;
		        gam1=chebev(c1,NUSE1,xx);
		        gam2=chebev(c2,NUSE2,xx);
		        gampl= gam2-xmu*gam1;
		        gammi= gam2+xmu*gam1;
		        ff=2.0/PI*fact*(gam1*System.Math.Cosh(e)+gam2*fact2*d);
		        e=System.Math.Exp(e);
		        p=e/(gampl*PI);
		        q=1.0/(e*PI*gammi);
		        pimu2=0.5*pimu;
                fact3 = (System.Math.Abs(pimu2) < EPS ? 1.0 : System.Math.Sin(pimu2) / pimu2);
		        r=PI*pimu2*fact3*fact3;
		        c=1.0;
		        d = -x2*x2;
		        sum=ff+r*q;
		        sum1=p;
		        for (i=1;i<=MAXIT;i++) 
                {
			        ff=(i*ff+p+q)/(i*i-xmu2);
			        c *= (d/i);
			        p /= (i-xmu);
			        q /= (i+xmu);
			        del=c*(ff+r*q);
			        sum += del;
			        del1=c*p-i*del;
			        sum1 += del1;
                    if (System.Math.Abs(del) < (1.0 + System.Math.Abs(sum)) * EPS) 
                        break;
		        }
		        if (i > MAXIT) 
                    throw new ApplicationException("bessy series failed to converge");
		        rymu = -sum;
		        ry1 = -sum1*xi2;
		        rymup=xmu*xi*rymu-ry1;
		        rjmu=w/(rymup-f*rymu);
	        } 
            else 
            {
		        a=0.25-xmu2;
		        p = -0.5*xi;
		        q=1.0;
		        br=2.0*x;
		        bi=2.0;
		        fact=a*xi/(p*p+q*q);
		        cr=br+q*fact;
		        ci=bi+p*fact;
		        den=br*br+bi*bi;
		        dr=br/den;
		        di = -bi/den;
		        dlr=cr*dr-ci*di;
		        dli=cr*di+ci*dr;
		        temp=p*dlr-q*dli;
		        q=p*dli+q*dlr;
		        p=temp;
		        for (i=1;i<MAXIT;i++) 
                {
			        a += 2*i;
			        bi += 2.0;
			        dr=a*dr+br;
			        di=a*di+bi;
                    if (System.Math.Abs(dr) + System.Math.Abs(di) < FPMIN) 
                        dr = FPMIN;
			        fact=a/(cr*cr+ci*ci);
			        cr=br+cr*fact;
			        ci=bi-ci*fact;
                    if (System.Math.Abs(cr) + System.Math.Abs(ci) < FPMIN) 
                        cr = FPMIN;
			        den=dr*dr+di*di;
			        dr /= den;
			        di /= -den;
			        dlr=cr*dr-ci*di;
			        dli=cr*di+ci*dr;
			        temp=p*dlr-q*dli;
			        q=p*dli+q*dlr;
			        p=temp;
                    if (System.Math.Abs(dlr - 1.0) + System.Math.Abs(dli) <= EPS) 
                        break;
		        }
		        if (i >= MAXIT) 
                    throw new ApplicationException("cf2 failed in besseljy");
		        gam=(p-f)/q;
		        rjmu=System.Math.Sqrt(w/((p-f)*gam+q));
		        rjmu=SIGN(rjmu,rjl);
		        rymu=rjmu*gam;
		        rymup=rymu*(p+q/gam);
		        ry1=xmu*xi*rymu-rymup;
	        }
	        fact=rjmu/rjl;
	        jo=rjl1*fact;
	        jpo=rjp1*fact;
	        for (i=1;i<=nl;i++) 
            {
		        rytemp=(xmu+i)*xi2*ry1-rymu;
		        rymu=ry1;
		        ry1=rytemp;
	        }
	        yo=rymu;
	        ypo=nu*xi*rymu-ry1;
	        xo = x;
	        nuo = nu;
        }
	    //void besselik(double nu, double x);

	    double jnu(double nu, double x) 
        {
		    if (nu != nuo || x != xo) 
                besseljy(nu, x);
		    return jo;
	    }
	    double ynu(double nu, double x) 
        {
		    if (nu != nuo || x != xo) 
                besseljy(nu,x);
		    return yo;
	    }
        //double inu(double nu, double x) 
        //{
        //    if (nu != nuo || x != xo) 
        //        besselik(nu,x);
        //    return io;
        //}
        //double knu(double nu, double x) 
        //{
        //    if (nu != nuo || x != xo) 
        //        besselik(nu,x);
        //    return ko;
        //}

        //void airy(double x);
        //double airy_ai(double x);
        //double airy_bi(double x);

	    void sphbes(int n, double x)
        {
	        double RTPIO2=1.253314137315500251;
	        double factor,order;
	        if (n < 0 || x <= 0.0) 
                throw new ArgumentException("bad arguments in sphbes");
           
            order = n + 0.5;
            besseljy(order, x);
            factor = RTPIO2 / System.Math.Sqrt(x);
            sphjo = factor * jo;
            sphyo = factor * yo;
            sphjpo = factor * jpo - sphjo / (2.0 * x);
            sphypo = factor * ypo - sphyo / (2.0 * x);
            sphno = n;            
        }

	    double chebev(double[] c, int m, double x) 
        {
		    double d=0.0, dd=0.0, sv;
		    int j;
		    for (j=m-1;j>0;j--) 
            {
			    sv=d;
			    d = 2.0*x*d - dd + c[j];
			    dd=sv;
		    }
		    return x*d-dd+0.5*c[0];
	    }

        static double SIGN(double a, double b)
        {
            return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);
        }
    }
}
